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Abstract 

In recent years many researchers have investigated the use of Markov random fields 
(MRFs) for computer vision. They can be applied for example in the output of the vi- 
sual processes to reconstruct surfaces from sparse and noisy depth data, or to integrate 
early vision processes to label physical discontinuities. Drawbacks of MRFs models 
have been the computational complexity of the implementation and the difficulty in 
estimating the parameters of the model. 

In this paper we derive deterministic approximations to MRFs models. One of the 
considered models is shown to give in a natural way the graduate non convexity (GNC) 
algorithm proposed by Blake and Zisserman. This model can be applied to smooth 
a field preserving its discontinuities. A new model is then proposed: it allows the 
gradient of the field to be enhanced at the discontinuities and smoothed elsewhere. All 
the theoretical results are obtained in the framework of the mean field theory, that is 
a well known statistical mechanics technique. A fast, parallel and iterative algorithm 
to solve the deterministic equations of the two models is presented, together with 
experiments on synthetic and real images. The algorithm is applied to the problem of 
surface reconstruction is in the case of sparse data. We also describe a fast algorithm 
that solves the problem of aligning the discontinuities of different visual models with 
intensity edges via integration. 
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1 Introduction 

In order to give a viewer information about a three dimensional scene many 
algorithms have been developed on several early vision processes, such as edge 
detection, stereopsis, motion, texture, and color . This information refers to 
properties of the scene as shape, distance, color, shade or motion, and it is 
usually noisy and sparse: more processing is then necessary to extract the 
relevant information and fill in sparse data. In recent years many researchers 
[9][12][6] [4] [5] have investigated the use of Markov Random Fields (MRFs) 
for early vision. MRFs models can generally be used for the reconstruction 
of a function starting from a set of noisy sparse data, such as intensity, 
stereo, or motion data. They have also been used to integrate early vision 
processes to label physical discontinuities. Two fields are usually required 
in the MRFs formulation of a problem: one represents the function that 
has to be reconstructed, and the other is associated to its discontinuities. 
The essence of the MRFs model is that the probability distribution of the 
configuration of the fields, given a set of data, is given as a Gibbs distribution. 
The model is then specified by an "energy function", that can be modeled to 
embody the a priori information about the system. In the standard approach 
an estimate of the field and its discontinuities is given by the configuration 
that maximizes the probability distribution, or equivalently that minimizes 
the energy function. Since the discontinuity field is a discrete valued field (it 
assumes only the values or 1) this becomes a combinatorial optimization 
problem, that can be solved by methods of the Monte Carlo type (simulated 
annealing[10], for example). 

The MRFs formulation is appealing because describes a system by local in- 
teractions and allows to capture many features of the system of interest by 
simply adding appropriate terms in the energy function. However, it has two 
main drawbacks: the amount of computer time needed for the implementa- 
tion and the difficulty in estimating the parameters that control the relative 
weight of the various terms of the energy function. 

In this paper we propose a deterministic approach to MRFs models. It 
consists in explicitly writing down a set of equations from which we can 
compute estimates of the mean values of the field / and the line process. 
We study two MRFs models, the second being an extension of the first, and 
we show that in each case the deterministic equations lead to a fast, parallel 
and iterative algorithm that can be used to obtain estimates of the fields 



everywhere. 

A natural framework for this approach is the equilibrium statistical me- 
chanics, since it deals with systems with many degrees of freedom described 
by Gibbs probability distributions. In principle, statistical mechanics allows 
us to derive the mean statistical values of the field (mean field equations) as 
explicit functions of the data and the parameters of the model: an algorithm 
implementing this would consist just of one step. Since the analytical com- 
putations required to obtain these explicit expressions are usually too hard, 
some approximations have to be made. We use the mean field approximation] 
a well known statistical mechanics tool, to obtain an approximated solution,' 
that is given in implicit form by a set of non linear equations. We call these 
equation deterministic to underline the deterministic character of the whole 
procedure. 

We concentrate our attention on the partition function Z, that is the sum of 
the probability distribution over all the possible field configurations, since it is 
known to contain all the information about the system. The idea underlying 
our approach is first to eliminate the line process degrees of freedom from Z: 
we will show that in so doing the effect of their interaction with the field can 
be simulated by a temperature dependent "effective potential" that depends 
only upon /. Its use is fundamental in the derivation of the deterministic 
equations, and gives useful insights on the role and the significance of the 
parameters. 

An advantage of such an approach is that the solution of the deterministic 
equations is faster than the Monte Carlo techniques, fully parallelizable and 
feasible of implementation on analog networks. The possibility of writing a 
set of equations is also useful for a better understanding of the nature of the 
solution and of the parameters of the model. 

We discuss two different MRFs models. The energy function of the first 
model has been already studied by several authors[2][ll][13][12]. It is inter- 
esting to notice that the GNC algorithm, proposed by Blake and Zisserman 
[2] in ad hoc manner, arises naturally in the framework of statistical mechan- 
ics. This estabilish a connection between MRFs and deterministic algorithm 
already used in vision. 

In the second model we define an energy function with an extra term, that 
establishes an interaction between the line process at neighborhood sites. 
This interaction can stimulate the creation of a line at a particular site if 
a line at a neighborhood site has been created. This term allows the data 



to be smoothed and, at the same time, the contrast at the discontinuities 
to be enhanced. Typical edge detection features like hysteresis, threshold 
and suprathreshold, which do not arise from the previous energy function, 
will arise naturally from the model We point out that this model is a gen- 
eralization of the previous one, that can be recovered by simply setting an 
appropriate parameter to zero. 

These two methods can be applied to dense data and, with small modi- 
fication, to sparse data as well. The problem of surface reconstruction (and 
image restoration) from sparse data is addressed and an algorithm to perform 
these tasks is obtained and implemented. We also outline an algorithm that 
solves the problem of aligning the discontinuities of different visual models 
with intensity edges via integration. 

The paper is organized in the following way: Chapter 2 presents an overview 
of MRFs in vision. Chapter 3 discusses the deterministic approximation of 
MRFs for the three energy functions mentioned above. Chapter 4 shows 
how to estimate the parameters of the models. In Chapter 5 some results are 
described. Chapter 6 discusses applications for sparse data and integration 
of visual modules with intensity. Chapter 7 concludes the paper. 

2 MRFs for smoothing fields and detecting 
discontinuities 

Here we briefly summarize how MRFs have been used in computer vision. A 
more extensive discussion is given in Geman and Geman [9], Marroquin[12], 
Chou[4], Gamble and Poggio[6], Gamble, Geiger, Poggio, Weinshall [5]. 
Consider the problem of approximating a surface given sparse and noisy 
depth data, on a regular 2D lattice of sites. We think the surface as a field 
(surface-field) defined in the regular lattice, such that the value of this field 
at each site of the lattice is given by the surface height at this site. The 
Markov property asserts that the probability of a certain value of the field at 
any given site in the lattice depends only upon neighboring sites. According 
to the Clifford-Hammersley theorem, the prior probability of a state of the 
field / has the Gibbs form: 

P(f) = Y, e - ms) (2-i) 



where / is the field, e.g. the surface-field, Z s is the partition function, U(f) = 
T,iEi(f) is an energy function that can be computed as the sum of local 
contributions from each lattice site i, and /? is a parameter that is called the 
inverse of the natural temperature of the field. If a sparse observation g for 
any given surface-field / is given and a model of the noise is available then 
one knows the conditional probability P(g\f). Bayes theorem then allows to 
write the posterior distribution: 

Pm - mmii s i e -^, 3 , . (2 . 2) 

As a simple example, when the surfaces (surface-fields) are expected to be 
smooth and the noise is Gaussian, the energy is given by 

E(f\a) = Ihi(/i - 9i? + «£(/, - m, (2.3) 

1 jcNi 

where 7,- = 1 or depending on whether data are available or not and N { is 
a set of sites in an arbitrary neighborhood of the site i. The maximum of 
the posterior distribution or other related estimates of the "true" data-field 
value can not be computed analytically, but sample distributions of the field 
with the probability distribution of (2.2) can be obtained using Monte Carlo 
techniques such as the Metropolis algorithm [14]. These algorithms sample 
the space of possible values of the surface-field according to the probability 
distribution P(f\g). 

One of the main attractions of MRFs models is that they can deal directly 
with discontinuities. Geman and Geman [9] introduced the idea of another 
field, the line process, located on the dual lattice, and representing explic- 
itly the presence or absence of discontinuities that break the smoothness 
assumption (2.2). The dual lattice is another lattice coupled with the data 
field lattice such that for each site of the data field lattice there are two sites 
of the dual lattice, one site corresponding to the vertical line and the other 
one to the horizontal line. The associated prior energy then becomes: 

*.■(/, = E (fi ~ mi ~ U 3 ) + E Veil*) (2.5) 

JeNi c 

where /,, is the element of the binary field / located between site ij. The term 
Vc(kj), where C is a clique defined by the neighborhood system of the line 



process (binary field), reflects the fact that certain configurations of the line 
process are more likely to occur than others. Depth discontinuities are usually 
continuous and non-intersecting, and rarely consist of isolated points. These 
properties of physical discontinuities can be enforced locally by defining an 
appropriate set of energy values V c {Uj) for different configurations of the 
line process ([9], [13], [6].) In our models the cliques will be simplified to the 
nearest neighbors. 

Two basic problems have arisen in using MRFs to solve vision problems. The 
first one is the amount of computer time used in the Metropolis algorithm 
or in simulated annealing [10]. The second problem is to how estimate the 
parameters of the energy function since they have been estimated in an ad 
hoc manner. 

We propose to approximate the solution of the problem formulated in the 
MRFs frame with its "average solution." The mean-field theory (which is 
explained in the next chapter) allows us to find deterministic equations ana- 
lytically for MRFs whose solution approximates the solution of the statistical 
problem. 

3 A deterministic approximation of MRFs 
3.1 The Line Process for two dimensions 

The line process was described in chapter 2 for the one dimensional case. 
Now we generalize it to two dimensions. In this case we define a horizontal 
line process h {j and a vertical line process v ir The line process h {j connects 
the site (ij) to the site (ij - 1), while Vij connects the site (ij) to the 
site (i - IJ). This is illustrated in figure 1. It is important to notice that 
the horizontal line process is determined by the gradient of the field in the 
vertical direction, and the vertical line process is determined by the gradient 
of the field in the horizontal direction. In this way we can reduce the two 
dimensional problem to two one dimensional problems, provided that the 
horizontal and vertical line processes do not interact. At the end we can add 
both line processes to get the edge map. This approach will be exploited in 
the next chapters. We point out that the fields f ih h {j and v {j are defined 
in the same lattice instead of having f {j in one lattice and h {j , v {j in a dual 
lattice. 
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Figure 1: The surface field f, the horizontal line process h and the vertical 
line process v are represented in two sites, (i,j) and (i+2, j+2), of the lattice. 
There are three fields defined in each site of the lattice as opposed to having 
one field in a lattice and the other two fields in a dual lattice. 



3.2 The Weak Membrane energy 

Or how to smooth the field but not at the discontinuities 
The Weak Membrane energy is given by 

E 1 (f,h,v) = E fg (f) + E Jl (f,h,v) + E,(h,v) (3.1) 

where 

E fAf) = J2(kj-9i,j) 2 (3.1a) 

E fl (f, h,v) = aJ2 Uj ~ Ai-i) 2 (l " h id ) + (fa - /,-_ lfi ) 2 (l - v itj )] (3.16) 



»»j 



Ei(h,v) = i'£ i ( h ij + Vij) (3.1c) 



*,3 



and oj and 7 are positive valued parameters. 

The first term, as in the previous case, enforces closeness to the data 
and the second one contains the interaction between the field and the line 
processes: if the horizontal or vertical gradient is very high at site (i,j) 
the corresponding line process will be very likely to be active (h { j = 1 or 
Vij = 1), to make energy decrease and signal a discontinuity. The third term 
takes into account the price we pay each time we create a discontinuity and 
is necessary to prevent the creation of discontinuities everywhere. 

3.3 Mean field theory and Weak Membrane 

We assume that there is uncertainty in the model and then the statistical 
framework can be used to incorporate the uncertainty as described by (2.1). 
Using mean field theory it can be shown (see Appendix A) that 

-j- 1 OF 

where F = -\\nZ is the free energy. Z is the partition function defined as 



in 

where J2{f} means the sum over all the possible configurations {/} of the 
system. 

The computation of the partition function is equivalent, in this case, to the 
evaluation of a multi-dimensional integral which can not be explicitly solved, 
due to the interaction between all the variables. Even if an exact solution 
can not be found, we can still obtain a good solution making use of the so 
called mean field approximation. The mean-field approximation is a general 
tool used in statistical mechanics that consists in substituting the interaction 
among the fields at different locations by the interaction of the field at each 
site with the mean field value at different locations (see appendix A) 
After this step the partition function factors into the product of single-site 
partition functions, that can usually be computed 

•■ {/> 

Now the partition function, and then the free energy are functions of the 
mean values of the fields, which are still unknown. These values, however, 
are usually related to the derivatives of the free energy with respect to some 
external field, as in (3.2). Substituting the free energy with its mean-field 
approximation, we obtain a set of implicit and self- consistent equations for 
the fields, usually called mean-field equations. In this sense the mean-field 
equations give a deterministic solution to the MRFs problem, since they 
relate explicitly the mean values of the field to the data. They could be 
implemented by a deterministic network. 
In the case of (3.1) the partition function becomes 

Z = £ c- /? S.-,>[(/.-.>-^.i) 2 +«(A? l /+A ? /)] j2 e -0Z i , j [ki.&, i +vi.jG!J (3>3) 

{/} {h,v} 

where G*. = 7 - aA?/, G^ = 7 - a^-\ Afc = fa - fa^ and Av. = 



3.3.1 Averaging Out the Line Process 

The contribution of the line process to the partition function can be exactly 
computed. Indeed the line process term in (3.3) is the partition function 
of two spin systems (h and v) in an external field (G h and G v ) with no 
interaction between neighboring sites. Then each spin contributes to the 
partition function independently from the others and its contribution is (1 + 
e~ pG i.j) for the horizontal field and a similar factor for the vertical one. The 
partition function can then be rewritten as 

Z = £ e^E^-K^-w^+^A^VM IJ(i + e -^) (1 + e -/^,) (3>4) 
if) ij 

The capability of computing the line process contribution to Z allows us to 
obtain a solution for the mean values of the fields h and v once a solution for 
the field is provided. As an intermediate step, after some algebra derived in 
appendix B we obtain 

hi=< TT7 s t > and %i =< i + \^i, > (3 - 5) 

where < ... > means the statistical mean value. In order to get a deterministic 
solution for the line process we need to do some approximations. We will 
assume then that we can neglect the statistical fluctuations of the field and 
so we replace the value of the field in (3.5) with its mean value. This is in 
essence the mean-field approximation. 

** = 1 + e nT-«Vu-Ti-i^) and ^ = ! + ^y-a^-ru-M ( 3 - 6 )' 

In the zero temperature limit (/? -» oo) (3.6) becomes the Heaviside function 
(1 or 0) and. the interpretation is simple: when the horizontal or vertical 
gradient (f itj - fa^ or fa - fa hj ) are larger than a threshold (^J) a 
vertical or horizontal discontinuity is created, since the price to smooth the 
function at that site is too high. This leads to a clear interpretation of the 
parameter 7, as it will be discussed in section 4.2. 

We notice that (3.6) can account for diagonal lines. We show with the 
example of a rectangle in a square lattice (see figure 2a). The threshold (J*) 




diagonal lines 
horizontal + vertical lines 



Figure 2: a) At the discontinuity, Aj;/ = (B - A) 2 and A? - 2 = 0. A 
chain of horizontal line processes are created for ( x < (B — A) 2 '), b) At the 
discontinuity, A?/ = (B - A) 2 and A?/ = (B - A) 2 . A chain of horizontal 
and vertical line processes are created for ( x < (B — A) 2 ). 



to create a discontinuity line is invariant under rotations of the rectangle with 
respect to the lattice (see figure 2b). The diagonal line is perceived from a 
chain of horizontal and vertical lines (like a stair case). The important result 
obtained from (3.6) is that although the total line created in figure 2b is \/2 
times bigger than in figure 2a the threshold to create a line has been kept 
the same. 

3.3.2 The Effective Potential 

We discuss how the interaction of the field / with itself has changed after 
the line process has been eliminated from the partition function. From (3.4) 
we notice that the partition function can be rewritten as 

Z = J2 e- 0{E fs>( f> > +E effW 

if} 
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Figure 3: The effective potential is shown as a function o/A*-. (For illustra- 
tion we set A£ = O.J a) For = 0.002. b) Zero temperature limit ((3 -+ <x>). 



where 



EefAf) = E«(4?/ + A"/) - I/„[(i + e-" G t)(l + e-Wj)] 

and Ef g (f) is given by (3.1a). 

This is the partition function of a system composed of one continuous valued 
field, whose energy is E fg + E eff . We interpret this result as the effect 
of the interaction of the line processes with the field /. This effect can 
be simulated by modifying appropriately the interaction of the field with 
itself, substituting the smoothing term in the energy function with a new 
temperature dependent potential. 

In figure 3 the effective potential is depicted for different temperatures. It 
simulates the effect of the line processes on the field /. Notice that the energy 
function is still the sum of local interactions between first neighbors. For the 
zero temperature limit one can see in figure 3 that the smoothing term is 
active only when the gradient is smaller than a threshold, proportional to 
the ratio between 7 and a. When the temperature is different from zero the 
border (threshold) of the smoothing region is no longer well defined due to 
thermal noise. It corresponds to an adaptive threshold that depends on the 
temperature. 
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3.3.3 A Deterministic Solution for / 

The introduction of the effective potential allows us to obtain a set of de- 
terministic equations for the field /. The mean field theory, as explained 
in section 3.1, can be applied to the partition function of equation (3.4) to 
obtain a set of mean-field equations depending on the temperature. In terms 
of statistical mechanics the data term plus the effective potential represent 
the free energy of the system. The mean field solutions are obtained by min- 
imizing the free energy. The set of deterministic equations can be written 
as 

d (E fg + E eff (f)) = 



and after some computation 

kj = 9i,j -_<x(fij_- Ai-i)(l ~ Vij) +_a(A;+i_- Aj)(l - v iJ+1 ) 

-<x(fi,j - fi- hj )(l - h itj ) + aifi+ij - fij)(l - h i+1<j ) (3.7) 

where hij and v iyj are given by formula 3.6. 

Equation (3.7) gives the field at site ij as the sum of data at the same 
site, plus an average of the field at its neighbor sites. This average takes in 
account the difference between the neighbors. The larger is the difference, 
the smaller is the contribution to the average. This is captured by the term 
(l-/ t)i ), where l itj is the line process. At the zero temperature limit (/? -> oo) 
the line process becomes 1 or and then only terms smaller than a threshold 
must be taken in account for the average. This interpretation helps us in 
understanding the role of the a and 7 parameters, as it will be discussed in 
chapter 4. Notice that the form of (3.7) is suitable for the application of a 
fast, parallel and iterative scheme of solution. 

3.3.4 The Effective Potential and the Graduate Non Convexity 
Algorithm 

We have to point out that this energy function has been studied by Blake 
and Zisserman [2], in the context of edge detection and surface interpolation. 
They do not derive the results from the MRFs formulation but they simply 
minimize the Weak Membrane energy function. From a statistical mechanics 
point of view the mean-field solution does not minimize the energy function, 
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but this becomes true in the zero temperature limit, so their approach must 
be recovered from the MRFs formulation in this limit. This is indeed the 
case, and it is easy to show that the effective potential becomes the Blake 
and Zisserman potential when f$ goes to infinity. In order to obtain the min- 
imum of the energy function E x Blake and Zisserman introduce the GNC 
(Graduate Non Convexity) algorithm, which is different from our determin- 
istic scheme, but can be embedded in the MRFs framework in a natural way. 
Let us review briefly the GNC algorithm. The main problem with the Weak 
Membrane Energy is that is not a convex function and a gradient descent 
method can not be applied to obtain the minimum because one could be 
trapped in a local minimum. In order to solve this problem Blake and Zis- 
serman introduce a family of energy functions E&\ depending continuously 
on a parameter p, pe[0,l], such that E™ is convex, E^ = E x and E^ are 
non convex for pe[0, 1). Gradient descent is successively applied to the energy 
function E^ for a prescribed decreasing sequence of values of p starting from 
p - 1, and this procedure is proved to converge for a class of given data. The 
construction of the family of energy functions JE^is ad hoc and uses piece- 
wise polynomials. In our framework, a family of energy functions with such 
properties is naturally given by E ef / T) where T is the temperature of the 
system. The GNC algorithm can then be interpreted as the tracking of the 
minimum of the energy function as the temperature is lowered to zero (like 
a deterministic annealing). In this way the approach of Blake and Zisserman 
can be viewed as a deterministic solution of the MRFs problem, even if it 
does not fully exploit the possibility of obtaining deterministic equations for 
the surface and discontinuity fields. The results obtained by applying this 
method to edge detection, for example, are good, and a pattern of meaningful 
discontinuities can usually be recovered [2]. However, sometimes the full set 
of discontinuities is not obtained: when the gradient of the image brightness 
is under the threshold, a discontinuity may not be detected, even though it 
would be necessary, for example, to close a contour. 

If interaction between lines (self interaction of the line-field) is introduced in 
the energy function this problem can be overcome, as we discuss next. 

Figure 4 summarizes the procedure used before to obtain deterministic 
algorithms from the statistical methods. 
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Sketch of the Method 



Construct Energy function 
(E depends on f and l) 



Build a statistical model 
Probability P = e -MS(f,l) / Z 




Mean field equations for f and 1 



More precisely , 



Compute Z and Mean field equations 



1. Sum P over 1 



2. Sum over f, but it is hard 



Obtain E e b ff (f) and 
in mean field ec 

Aditional feature: Deterministic annealing (increase b) 



So, find f that minimise E^f ) 
obtain mean field equation for 1 (mean field approximation) 



Figure 4: Overview of the method 
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3.4 Improving the weak membrane model 

Or how to smooth the field while enhancing the differences at the discontinu- 
ities 

So far we have not exploited an important physical constraint of images, 
namely the smoothness of the discontinuity field. Isolated discontinuities are 
very unlikely to occur and, on the contrary, the presence of a discontinuity at 
a site makes more likely the presence of a discontinuity at a neighboring site. 
This smoothness constraint on the discontinuity field can be incorporated in 
the model by simply adding a new term to the energy function. Thus, the 
total energy becomes 

E 2 = Ei -f Eu 
where E x is given by (3.1) and we define the new term 

En = -tlY^KlKj-i + V iJ v i-hj] 

and e is a new parameter whose exact meaning and estimation will be ex- 
plained in the next section. To make more evident the meaning of the new 
term we notice that 



E, + E u = £ 7(1 - 2e)[/i? J + v^] + e 7 p,. ,- - ft,, ,_ x ) 2 + (»„ - «,_,„•)»] (3.8) 

where E t is given by (3.1c). From (3.8) it is clear that e is related to the degree 
of smoothness of the discontinuity field and so must be a positive number. 
In order to keep positive the price we pay for creating a discontinuity and 
to prevent the line processes from being active everywhere, t has to be less 
than 1. 

We notice here that this model is a simplified version of other possible po- 
tentials. In particular, the neighbor size considered here is at most of two 
pixels where Gamble k Poggio (1987) for example, have discussed more so- 
phisticated cliques composed of larger neighbors. 
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3.4.1 Averaging Out the Line Process 

Again, as in the Weak Membrane energy case we calculate the contribution 
of the line process in the partition function. The partition function can be 
written as 

x E { »,. } e-^.. J ^^^- ii -^^) + « i .(^-e,^i^)] (3 . 9) 

where the Gs have been defined in section 3.3. We rewrite the line process 
term by substituting h itj (G^j - erf hi j^) by h itj (G^ - C7 i*ii=l±^i±i). Notice 
that this change is just a rearrangement of the terms on the sum over i,j and 
therefore does not change the value of Z. We also notice that the line process 
term of (3.9) is the partition function of many one- dimensional Ising spin 
models in an external field. For a constant external field an exact solution 
can be obtained by the method of the transfer matrix (see Parisi [15]). This 
method can be still applied to our case (not constant external field), giving 
a compact expression for the partition function. However, the expression we 
obtain is highly non local for the field / and not useful for fast computations. 
Because of it we_ apply the mean field approximation and replace A f j_i by 
hi,j-u hij+i by h itj+1 , u t -__ lti by t) t _i,j, and v i+1 j by v i+1J in (3.9). Notice 
that at this step the mean field approximations refer to the line process and 
not to the field / (that is kept "frozen"). The sum over the line process 
configurations (3.9) can then be computed, giving 

Z mf-hv = £ c "/'E.-,iK/*.i-W.i) 2 +o-(Aj i 2 +A? >i a )] x 

The partition function is now similar to the previous one (3.4) and the zero- 
temperature mean-field equations for the line process field can then be de- 
rived in the same approximation: 

hi = MdAi - A,-:) 2 - 7 + c 7 ^t±±hi±l) and 



2 



«« = "/»(«(/« - fi-i,i) 2 - 7 + rt f>, - 1J t ,><+1J ) 



2 



(3.10) 



16 



where <Tp(x) = ^^ is the sigmoid function. Notice that now equations 
(3.10) form a set of non linear equations and the solution for the line process 
is not as simple as before. 

3.4.2 Hysteresis, Suprathreshold and Threshold 

The solution obtained in (3.10) for the line process can deal with the prob- 
lem of "streaking". Streaking is the breaking up of an edge contour caused 
by fluctuations above and below a threshold along a contour [3]. This is a 
common problem with thresholded detectors, in particular with the scheme 
implemented by the Weak Membrane energy. For the third energy the thresh- 
olding is done with hysteresis. This is the way that Canny's detector works 
[3]. The hysteresis phenomenon is evident in (3.10). With the creation of a 
horizontal line at a neighbor site, say {ij - 1} {hij^ = 1), the energy neces- 
sary for creating a line at a site {z, j] decreases by f. On the other hand, if 
two lines are created at {ij - 1} and {ij + 1} then the energy decreases by 
67. A low threshold (threshold) and a high threshold (suprathreshold) arise 
naturally. The suprathreshold for creating a line is given by ^ < (Aj 1 ) 2 , in 
this case a line is al ways created. In the same way the lowest threshold is 
given by J(l - e) < (A^) 2 < J, in this case a horizontal line will be created 
at site {i,j} if a horizontal line has been created in the sites {i,j - 1} and 

0\j + i}. 

At higher temperatures the threshold and suprathreshold are not so well - 
defined, and what we have are adaptive thresholds that depend on the values 
of the field. These points will become clearer when we discuss the effective 
potential. 

We point out that some hysteresis may occur on the Weak Membrane en- 
ergy. However it will happen just on special cases and without much of 
control. A consistent framework to enhance and smooth images requires at 
least two thresholds. The Weak Membrane energy is a one-threshold model 
and therefore can not fully exhibit hysteresis. 

3.4.3 The Effective Potential 

By analogy to section 3.3, the effective potential is given by 
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EeSl(f) = Ei,i{«(A?/ + Av/) 

-I/n[(l + e-^tr ^ 1 - 1 ?-* 1 ))^ + e -/»(^- g > --'-1-t.rf ))] | (3 n) 

An exact computation by means of the transfer matrix method shows that 
this potential is highly non local and its practical use is limited. What we 
will do next is to use approximation techniques to obtain a potential that is 
local , up to the first neighbors. 

Without losing any generality we now dismiss the contribution of the vertical 
process since it is independent of the horizontal process contribution. In 
order to obtain an expression for E eff in terms of the gradient field A^ we 
investigate (3.10). Since the discontinuity field is smooth we first approximate 
the term m-i+*m-m by h {j and (3.10) at the zero temperature becomes 

hi = 0(<*(fij - fij-i) 2 - 7 + elks) (3.12) 



Ki = < 



The solution of (3.12) depends on the value of A^ = (fij-fij-i) 2 as follows 

' if (A£) 2 < 5(1 - e) 

Oorl if J(i _ e ) < (AJT)2 < x 

These results are illustrated in figure 5. We notice that in the region 2 (1 - 
e ) < i^ij) 2 < J both solutions are possible , thus reflecting a weakness in 
the assumption we made. A possible way to obtain a unique solution is to 
consider an average of the two solutions such that the higher the gradient of 
/, the higher is the value of h (more likely for creating a line). We interpret 
the values of h between 1 and as the likelihood of creating a line. We ap- 
proximate the solution of the mean field equation for h itj at any temperature 
(/?) by the analytical expression 



hi = ^(A*.) 2 - ^/n(l + e-<*—<^> 2 >) + JL/»(i + e -*) (3 . 13) 



18 




■>■ ij-1 



*• "ij.j 



Figure 5: The three figures illustrate the solution of (3.12) according to dif- 
ferent values of the gradient field f (A£) a) The line process is on (h^ = 1) 
b) Two solution (three solution for temperature different then zero) c) Shows 
the solution of no line-process (hij = 0) 






Figure 6: h iyj as a function of the gradient field (A*-). Two solutions of 

(3.13). a) /3 = 1 b) (3 = 0.002 



19 





-at -10 



to 20 



Figure 7: The third energy potential for /3 = 0.002 and different values for the 
parameter t. a) For e = 0.9 b) For e = the potential becomes the effective 
potential of the Weak Membrane energy. 

that has the desired property and should be a reasonable guess for the de- 
pendence of hij on the gradient of the field /. 

Figure 6 shows two solutions of (3.13) for different temperatures. Substitut- 
ing the value of h {j given by (3.13) in (3.11) we then obtain the third energy 
potential, that is plotted in Figure 7 for some values of /? and e. 
The effect of the third energy potential on the mean field / can be understood 
by analyzing the "force" (derivative) of the potential (see Figure 8). The 
third energy effect can be described as follows: 

1. Smooth the field for small values of the gradient ((Aj 1 ) 2 < (A ) 2 ). 
(Positive derivative of the potential.) 

2. Enhance the discontinuity for (A ) 2 < (A*-) 2 < J, in other words, 
forcing the value of the field at a site to be pulled away from the value 
of the field at the considered neighbors. (Negative derivative of the 
potential.) 



3. Neither smoothing nor enhancing the field for (A- 1 -) 2 > 
derivative of the potential is zero.) 



(The 



We point out that the third energy potential derived above is a local approx- 
imation of the mean field approximation over the discontinuity field. 
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Figure 8: The third energy force for (3= 0.002 and different values of the 
parameter e a) For t=0.9 . b) For e = the force is the same as for the Weak 
Membrane energy. 

3.4.4 A deterministic solution for / 

Here again the same arguments of section 3.3.3 apply to obtain a determin- 
istic solution for the field /. After some computations we then obtain: 



J*J — 



g^-a/^Ul - Vij) + aAl i+1 (l - v iJ+1 ) 
-aAj^l - hij) + aAg^l - h i+ ij) 
+ ot^jhiMh - *(Afo)3 

- acAk +1 ft tW <rfl(7 - a(A»,- +1 ) 2 ) 
+ aeAl-Vija^ - afAv,.)^) 

- acAv i+1 i; <ii+1 a- / j(7 - a(Ay J+1 ) 2 ) 



(3.14) 



4 Parameters 

The parameters a, 7, e and (3 must be estimated in order to develop an 
algorithm that smoothes, enhances and finds the discontinuities of the given 
data-field. 

4.1 The parameter a 

The parameter a controls the balance between the "trust" in the data and 
the smoothing term. The noisier are the data the less you want to "trust" it 
so a is larger, the less noisy are the data the more you "trust" it so a should 
be smaller. To estimate a various mathematical methods are available. The 
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generalized cross validation method introduced by Wahba [17] and the stan- 
dard regularization method described by Tikhonov [16, 1] were studied by 
Geiger and Poggio [7] to estimate a. The standard regularization method 
gives [7] an equation relating a with the signal-to-noise ratio that for the 
continuous case is 

f°° ["> S(u,v)(u; 2 + v 2 ) 2 J too ,«> N(uv)(lo 2 + v 2 ) , , , x 

where 5(u>,i/) is the spectral density (power spectrum) of the signal f(x,y) 
and N(w, v) is the spectral density of the noise v(x, y), assuming that v(x, y) 
is stationary. 

When an estimation of the noise is not available, it is necessary to use 
the generalized cross validation method (GCV). GCV states that the opti- 
mal value of a can be obtained by minimizing the functional (here in one 
dimension) 

\/(^\ — V^ [jn,a\U) — 9i\ 2/ \ i A rt \ 

v(a)_ »£ (i -«»(«))* Wt(a) (4 - 2) 

where /„,«(*,-) is the smoothed solution given by (3.4), 

u*(a) = (1 - a kk (a))/(l - i £ «;;(<*)) 



n j=o 



and a kk (a) = ^(/ n ,a)(^). 

For the purpose of smoothing images both methods give about the same 
values of a [7]. Better results should be obtained if the estimation of a is 
done locally (for a small neighborhood). 

4.2 The parameter 7 

From (3.6) one can see that y^J is the threshold for creating a line in the 
Weak Membrane energy. Let's call f the value y^J. From the expression of 
the effective potential we notice that if the gradient Afy is above £ there is 
no smoothing and if the gradient is below f then smoothing is applied. The 
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parameter £ defines the resolution of the system and we explain this with 
two examples. 

In the first example suppose we are working with the stereo module, so 
the data field is a depth-field. In this case £ is the threshold for the changes 
in depth to be called a depth discontinuity. This value is determined ac- 
cording to the resolution of the stereo system available and /or to the desired 
resolution that one is interested. 

In the second example the field is the intensity data and the parameter £ 
represents the threshold for detecting edges. This value is somehow arbitrary, 
and probably context dependent. A variation of £ = 16 units to f = 32 
units for an 8-bit array image is likely to be in agreement with the human 
threshold. The exact value of £ depends on the attention of the observer 
and/or the sensitivity of the system. 

For the second energy the absolute threshold to create an edge is also 
given by £. However when there is support from neighbor edges this threshold 
can be lowered to £ x (y^l - e) (see (3.10)). This suggests setting 7 so to 
guarantee that the highest noise gap is smoothed (we are assuming that noise 
gaps are isolated features). In this case the e parameter will be set to assure 
that at the edges the images will not be smoothed but perhaps enhanced. In 
this case the value of £ = 30 may be desired for an 8-bit array. 

4.3 The parameter e 

The parameter e allows the energy to be more general by controlling the 
amount of propagation of the line. So, once a line is created, the price to pay 
for creating another line next to it will be lowered by the amount of 7c In 
other words, from (3.11) one can see that the difference in the energy for when 
a line has been created and when no line has been created at pixel (i — 1, j), 
is given by 7c. This is what characterizes the threshold and subthreshold 
or the hysteresis phenomena [3]. The threshold is given by \/^f^ and the 
suprathreshold by ^/j. The parameter e varies from to 1.0 and when is zero 
reduces the third energy to the Weak Membrane energy. When e = 1.0 lines 
are created everywhere, since once a line is created there is no cost in creating 
another one and then it propagates indefinitely (to the image limits). How 
much does one want to propagate a line? How much should the difference 
between the threshold and the suprathreshold be? Which exact value of e 
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should be chosen? According to the discussion above about 7 we conclude 
that e should be chosen (assuming we already have a and 7) to guarantee 
that y 7 ' 1 " 6 ' is below the desired edge threshold. Again it depends on what 
one wants to do with the data. For detecting edges of objects in a scene one 
wants to have high suprathreshold and threshold (usually object boundaries 
exhibit high gradients) and e large (bigger than 0.5) so that all the object 
boundaries are detected, included the exceptional boundary pixels with a 
smaller gradient. 

4.4 The parameter (3 

The parameter (3 controls the uncertainty of the model. The smaller is /? 
the more inaccurate is the model. This suggests that for solving the mean 
field equations a rough solution can be obtained for a small value of j3 (high 
uncertainty) and thereafter we can increase /? (small uncertainty) to obtain 
more accurate solutions. This can be called deterministic annealing. 
We also notice that /3 multiply the first term on (3.1) to give 

which suggests that the inverse of /? gives the standard deviation of the noise. 
We keep in mind that /? also multiplies the other parameters of the model 
and therefore by estimating f3 by the amount of noise the estimation of the 
other parameters have to be scaled by this factor. 

5 Results 

5.1 Implementation 

For the implementation the zero temperature limit equations have provided 

results as good as the deterministic annealing with a faster computational 

time. We do not have proof of convergence but only suggestive experimental 

results. 

In order to find the mean field solution for the third energy, we solved (3.14) 

together with (3.10) in a coupled and iterative way. More precisely: 

For simplicity, let's write (3.14) as 
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JhJ — * \9i,ji Ji,j-> fi-l,j-> fi,j-li fi+l,j j fi,j+li ^t\j» ^i-l,ji v i,ji v i,j-l) 

and (3.10) as 

hj = H (fi,h fij-u hij-u hi,j+i) v itj = V(f itj , fi-ij, f t _i,j, v i+ i fj ) 

The algorithm solves these equations iteratively by updating (for the parallel 
case) in the following way 

f? 1 = %-"% - ^(ftj.&. •..,&«,*?,,... ,t£_0] (5.1) 

and 

KJ 1 = mft^fT-u, *&-!,*?.,•+,) and v#* = V(fTj, JTj-u Cj, v^j) 

where u controls the rate of change and the index n indicates the step of the 
iterative process. 

For a serial implementation we first update the even sites (like the white 
squares in a chess board) and then the odd sites (like the black squares). 
Typically the algorithm has converged in 10 iterations which takes about 1 
minute for images of 64 x 64 pixels on a Symbolics 3600. 

5.2 Synthetic Images 

We first tested the algorithm on a synthetic step edge image with noise added. 
This is a good test-image since locally many edges on real images are of this 
type. 

The step edge image is an 8-bit array of 64 x 64 pixels with a step intensity 
of 140 units (see Figure 9a). White noise, with standard deviation 30, is 
added to the step edge (see Figure 9b). We apply the algorithm described by 
(5.1) to reconstruct the original image. First we use e = 0, in order to have 
the Weak Membrane energy. We set a = 4 (typical smoothing parameter) 
and 7 = 24000 (so that we do not smooth the edge). The result after 20 
iterations is shown in Figure 9c. We notice that the step edge starts to be 
smoothed before all the noise has been smoothed. In Figure 9d we set e to be 
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0.9 and 7 = 45000, and the result is more striking. Because we set the ratio 
J to be higher than that, all the noise has been smoothed away. However, at 
the edge there is support from the neighbor edges to lower the threshold to a 
factor (1 - e) = 0.1. Therefore the noise is suppressed and the edges are not 
smoothed but on the contrary are enhanced. We conclude that for a noisy 
step edge, the third energy is able of retrieving and enhancing the edge (see 
figure 9 Id) while the Weak Membrane energy retrieves the edge with worse 
performance (see figure 9c). 

5.3 Real Images 

When we apply the second energy algorithm to a real still life image the 
result is an enhancement of specular edges, shadow edges and some other 
contours while smoothing out the noise (see Figure 10). This result is ob- 
tained consistentently in all the images we used and also with other fields 
besides intensity, e.g. "color" images. In our case "Color" images are 8-bit 
images representing the ratio between the red array and the sum of red and 
green arrays. 

6 Surface reconstruction and alignment via 
integration 

6.1 Sparse Data and Surface Reconstruction 

Until now we dealt with dense data defined on allthe sites of the lattice. 
This is not always the case, especially if the field is the output of a stereo or 
motion module. In that case data are given on a subset of the lattice sites 
and the data-field interaction term (Ej g ) in the energy function has to be 
modified in the following way 

hi hj 

Here 7^ is a flag that is one if a datum is given at site (i,j) and 

otherwise. This slight modification has no effect on the theoretical results, 

and only some changes take place in the deterministic equations for the 

surface field: the term enforcing closeness to data disappears when a datum 
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p d. 




Figure 9: a) Step edge with 140 grey value units for the step, lb) White 
noise with standard deviation 30 grey value units has been added, c) The 
noisy image after 20 iterations for a = 4, 7 = 24000, e = 0. d) The noisy 
images after 20 iteration for a = 4, 7 = 45000, e = 0.9 
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Figure 10: a) Still life image 128 x 128 pixels, b) The image smoothed with 
e = 0.5 ; 7 = 1400 and a = 4 for 9 iterations c) The line process field (without 
thinning). 

is not given at a site. We rewrite here the deterministic solution for the field 
/ in the case of the Weak Membrane Energy, since the third energy equations 
are modified in the same way. In the case of sparse data (3.7) becomes then 

hjH,J = 9i,j7i,j ~ <x(fij - fi,j-i)(l - Vij) + a{fij +1 - /^)(1 - vij+i) 
-<*{hj ~ /.-i,j)(l - Ki) + <*(fi+i,j ~ fi,j){l - h i+ u) 

To apply this algorithm one needs first to fill in the data. We chose to fill 
in the data by averaging next neighbors and applying the algorithm at the 
same time. So at each step of the algorithm each lattice site is visited: if 
there is no field value the average neighbor value is taken; otherwise we apply 
the above algorithm. We notice that no action is taken if at a particular site 
there is no field value and no neighbor field value. 

From one face image we produced sparse data by randomly suppressing 70 
% of the data (see Figure 11). We then applied the third energy algorithm 
to sparse data. The parameters were kept the same as the other real image. 
The reconstruction from sparse data can be applied to depth data in which 
case it is usually called surface reconstruction. We used a stereo algorithm 
based on zero crossings to obtain depth data from the image shown below. 
Then we applied the algorithm to reconstruct the depth surface. The param- 
eters were different according to the criteria for depth discontinuity. 
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Figure 11: a) A face image of 8-bits and 128 X 128 pixels, b) Randomly 
chosen 70 % of the original image. For display the other 30% are filled with 
white dots, c) The algorithm described above is applied to smooth and fill in 
at the same time with e = 0.9, 7 = 1400 and a = 4 for 70 iterations. 
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Figure 12: a,) The left image (8-bits and 128 x 128 pixels), b) The right 
image of the same scene c) Sparse depth data obtained by the stereo algorithm 
based on zero crossing. Intensity represents depth values, d) The algorithm 
is applied with e = 0.0, 7 = 400 and a = 4 to reconstruct the surface, e) 
Canny edges for figure a), f) The alignment algorithm is applied with the 
same parameters as d) and 8 = 384. Depth discontinuities aligned with the 
intensity edges are obtained. 
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6.2 Alignment of visual modules with intensity edges 

The integration of different visual modules to improve the detection of the 
discontinuities can also be addressed in this scheme. As suggested by Gamble 

6 Poggio [6] , we can add the term £ (v,-j + /i,j)(l -e^) to the Weak Membrane 
Energy or the second energy. Here e tJ is an external field, for example the 
edge map that is coupled with the stereo field. For implementation purposes 
the only consequence of adding this term is the change of the global parameter 

7 into the local parameter 7^ = y - 6(1 - e tj ). R. Thau implemented 
this scheme in the Connection Machine using Canny's intensity edges for 
e,j. We first took a pair of images (Figure 12) and used a stereo algorithm 
to find depth at the zero crossing positions. We then used the algorithm 
to reconstruct depth everywhere and detect depth discontinuities. As we 
expected depth discontinuities are aligned with the intensity edges. 

7 Conclusion 

We have used statistical mechanics tools to derive deterministic approxi- 
mations of Markov random fields models. In particular we have studied an 
energy model that is suitable for image reconstruction or any field reconstruc- 
tion. The model has been developed to include the following characteristics: 

• the surface field is smoothed when its gradient is not too high, 

• contrast will be enhanced where a discontinuity occurs (if it is not too 
large already), 

• the discontinuity field is likely to be smooth (isolated discontinuities 
are inhibited), 



• hysteresis and adaptive multiple threshold arise naturally from the 
model. 

• three parameters are needed to specify the model, 

• when one of the parameters (e) is set to zero the weak membrane energy 
is recovered. 

• An understanding of the role of the parameters is possible, 
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• The model can deal with sparse data and alignment of the discontinu- 
ities of different modules with the intensity edges. 

We have shown that the deterministic algorithm of GNC can be regarded 
as an approximation of the gradient descent method with a deterministic 
annealing schedule to solve the mean field equations. This suggests a unified 
framework to connect different methods used on image segmentation, restora- 
tion and surface reconstruction. We will show in another paper [8] that several 
deterministic algorithms for image segmentation and reconstruction are ap- 
proximations of two methods to solve the mean field equations: the gradient 
descent method discussed in this paper and the parameter space method 
discussed in [8]. We derived a deterministic solution for the mean values of 
the surface and discontinuity fields, consisting of a system of coupled non- 
linear equations. An algorithm has been implemented to obtain a solution 
for this system: it is fully parallelizable, iterative and recursive, allowing effi- 
cient computation. It would be interesting to analyze other approximations 
or extensions of this model. For example, a model that includes interac- 
tion between the horizontal and vertical line processes could be developed to 
inhibit self- intersections of the discontinuity field. A term like hij(l — Vij) 
may be sufficient. An analysis of convergence of the algorithm would also be 
important. 
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Appendix A 

Given an energy function of the form 

£(/) = £{(/<- ft) 2 + [£V(/,-,/ j )]} 

* jtNi 

where JV t - is the set of neighbor sites of i, V(fafa) is an interaction potential 
and g { is an external field (data), the partition function is given by 

f 
where the J2f means a sum over all the possible configurations of the system, 
that is a multidimensional integral over all the variables /,-, in this case. The 
mean field equation for fa is given by 

L f 

From this definition, and the definition of the partition function Z , the 
following equality is derived: 

zWk = ' 2fi(fk ' 9h) {A2) 

Defining as usual the free energy FasF= -jlnZ we can now obtain, from 
(A2) 

j IdF 

In the case of energy E the mean-field approximation consists in replacing 
the fields at the neighbor sites with their statistical mean value fa. Therefore 
the mean-field approximation E mf to the energy function becomes 

£■"'(/) = £{(/*- so* + ie ^ (/<>/;•)] 

» jeNi 

This is a good approximation of the original energy when the fluctuations of 
the field / are small, which is the case in our experiments. 
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Appendix B 

In this appendix we will derive an expression for the mean value of the h 
and v fields of a system described by the Weak Membrane Energy 

E,(f, h,v) = E„(f) + E t ,(f,h,v) + E,(h,v). 
Noticing that 

£,(/, h,v) = E 1 (f) + YllhjGli + VijGtj] 

where the energy function E Y has been defined in section 1, the partition 
function can be written as 

{/} {h,v} 

Let us consider just the field h. By definition the mean value h^ is 



* if) ih,v} 

Then we can rewrite (Bl) in the following way: 



(Bl) 



h^-^e-^^-ZUf) (B2) 



where 



Z hv (f) = E e-'Za^^WJ. 

{h,v} 

Zhv represents the contribution of the line process to the partition function 
of the whole system for each fixed configuration of the field / and can be 
seen as the partition function of the line process system when the field / is 
kept "frozen" . We notice that the line process system is simply a classical, 
non-interactive spin system in an external field. Since there is not interaction 
between the spin variables, Zh v {f) can be easily computed, giving 

i,3 
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We can now define the free energy of this spin system by the usual relation 

Z hv (f) = e-<^«<» 
and substituting in (B2) we obtain the following equation 

Before proceeding further we now notice that, having defined F hv (f), the 
partition function of the whole system can now be rewritten as 

2 = r e -^(/)+F^(/)]_ 
{/} 
This can be regarded as the partition function of a system, composed just by 
the field /, whose energy function is E x (f) + F hv (f). Now, indicating with 
< ... > the statistical mean value over all the possible configurations of this 
system, (B3) can be rewritten now in the following way: 

i.. = <2M > 
,J < dGfj > 

and then 

^ =< 7— W > • < B4 ) 

1 + er *>j 
A similar expression holds for Vij. The meaning of these equations is the 
following: if the field / is kept "frozen" the mean value of the horizontal 
line process is simply d ^i f ^ , but to obtain the true mean value hi -• we have 

to average d ^ h ^ over all the possible configurations of the field /, whose 
probability distribution is given by 

e [Ei{f)+F hv (f)] 

where the term Fh v {f) is due to the interaction of the line process with the 
field /. 
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